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The aim of the current work is a detailed time-series analysis of the data from Dark Matter 
direct detection experiments as well as related datasets. We examine recent claims that the cosmic 
ray muon flux can be responsible for generating the modulation signals seen by DAMA and, more 
recently, by the CoGeNT collaboration. We find no evidence for such a strong correlation and 
show that the two phenomena differ in their power spectrum, phase, and possibly in amplitude. In 
addition, we investigate in more detail, the time dependence of Dark Matter scattering. Since the 
signal is periodic with period of a year (due to the Earth's motion around the Sun), the presence 
of higher harmonics can be expected. We show that the higher harmonics generically have similar 
phase to the annual modulation and the biannual mode in particular could provide another handle 
in searching for Dark Matter in the laboratory. 



I. INTRODUCTION 

Whereas the gravitational nature of Dark Matter (DM) 
is a crucial ingredient for the success of the standard cos- 
mological model, its non-gravitational character proves 
elusive and remains essentially unknown to date. The 
possible connection with the electroweak scale and the as- 
sociated expectation of an interaction strength not much 
smaller than Gp and a DM mass not much lighter than 
0(10 GeV) has motivated the construction of experi- 
ments looking for the direct scattering of DM against 
nuclei in the lab [1] . The experimental efforts in the field 
are proceeding in full force and at least another decade 
of progress is expected. 

With the main observable being the nuclear recoil spec- 
trum the information content is rather limited. More- 
over, as was realized over the past decade, the recoil 
spectrum resulting from DM collisions is model depen- 
dent and may not follow the simple exponential rise to- 
wards low energies expected from elastic scattering ;2;- 
|6]. On the other hand, one of the most robust predic- 
tions of the cold DM model is that the relative velocity 
of the Earth against the DM halo should vary with the 
Earth revolution around the Sun [3 H]. The expected 
recoil rate is therefore a general periodic function with a 
fundamental period of one year with a particular phase. 
Despite the strong motivations to look for such modu- 
lations, it has so far been achieved by only two experi- 
ments, DAMA [HI Ho] and CoGeNT [U]. The difficulty 
is of course in maintaining the stability of the appara- 
tus and the rejection of background over a long period of 
time. Indeed, the two experiments just mentioned have 
experienced considerably more background events than 
some of the other extremely clean experiments such as 
CDMS-II [Hin] and XENONIOO T^. Moreover, the 
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search for annual modulations is not without difficulties 
since many more mundane phenomena arc known to ex- 
hibit such modulations. Therefore, two of the main goals 
of this paper are 1) to show that the time series of the re- 
ported DAMA and CoGeNT signals can be used to make 
firm statements about background hypotheses proposed 
in the literature and 2) to establish further observables 
based on the time distribution of the energy spectrum. 

The only direct detection experiment which claims de- 
tection of a firm signal from dark matter is DAMA, situ- 
ated in the underground LNGS laboratory at Gran Sasso, 
Italy. The DAMA dataset consists of two main periods, 
DAMA/Nal (Dec 1995 - July 2002) and DAMA/LIBRA 
(Sept 2003 - Sept 2009), amounting to a cumulative ex- 
posure of 1.17 tonxyears. The residual rate shown by 
the collaboration exhibits a very clear annual modula- 
tion compatible with what is expected from DM models 
where the Earth's motion around the Sun results in a 
modulation peaking on approximately June 2"^*. The col- 
laboration has also released the modulation amplitude as 
a function of recoil energy, which seems to exhibit a pe- 
culiar form possibly more consistent with inelastic scat- 
tering [2HS]- It should be kept in mind, however, that 
this energy spectrum is obtained from the full data set 
by assuming the periodic functional form. Sadly, there 
has been no release of the full data set to date. It is 
particularly unfortunate because, as we will see below, 
the procedure carried by the collaboration to obtain the 
power spectrum of the signal and other parameters makes 
it difficult to compare to both background and signal hy- 
potheses. 

Given the seasonal nature of the signal it is en- 
tirely conceivable that environmental effects induce back- 
grounds with an annual modulation pattern just like the 
one seen by DAMA. As such, it is clear that such contam- 
ination may depart significantly from a sinusoidal dis- 
tribution in time. It is therefore important — whenever 
enough data on a putative background inducing process 
is available — to assess its viability based on a full time- 
series analysis. For example, in this work we employ 
Pearson's coefficient of correlation as a test statistic when 
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FIG. 1. Percent annual residuals of the LVD measured muon 
fliuc when binned in accordance with the DAM A/LIBRA 
runs 1-5. The latter residuals are shown for the 2-4 keV bins 
assuming a baseline count rate of s = 1.15 cpd/kg/keV. 



quantifying the compatibility of two data sets. This al- 
lows for a quantitative comparison between two data sets 
which is independent of any assumption about the func- 
tional form of the time variations. 

Fluxes of neutrons probably constitute the most dan- 
gerous source of background as they lead to nuclear re- 
coils which can mimic DM- nucleus interactions [THj. In- 
deed, it is known from ICARUS measurements at LNGS 
[16j that the ambient neutron flux generated in the sur- 
rounding rock shows some variation on the timescale of 
a year. However, with a total of five data points we find 
that it is not possible to make any statistically signifi- 
cant statement regarding the temporal compatibility of 
this candidate background with the DAMA signal. We 
therefore choose not to elaborate on this issue any fur- 
ther. 

In contrast to rock-generated neutrons, a wealth of 
data is available on another guaranteed source of neu- 
trons: cosmic ray muons which can penetrate deep un- 
derground and induce spallation reactions in the detector 
and nearby. It is also possible for these muons to de- 
posit their energy directly into the detector crystals [l7j . 
Measurements of the muon flux underground have a 
long history and its seasonal variation is firmly estab- 
lished. There are published measurements at LNGS from 
MACRO pm^, LVD [19] and most recently from Borex- 
ino [20] • For DAMA, the relevant measurement is the one 
from LVD since it was taken at the same underground 
lab and it overlaps in time with the DAM A/LIBRA runs 
1-5. Figure [T] shows the percent residuals of the muon 
flux when binned in concordance with DAMA with the 
annual mean count rate subtracted. Also shown are the 
residuals in the 2-4 keV bin of the DAM A/LIBRA data 
assuming a baseline rate of s = 1.15 cpd/kg/keV. The 
seeming similarity in time and amplitude is tantalizing 
and it has lead various authors to put forward the hy- 
pothesis that both signals are in fact measurements of 
the very same cosmic ray phenomenon |151 1171 121] . It 
will be one of the central points of the paper to critically 
examine these claims, finding that they have difficulty 
explaining the observed modulations. 

The issue of annual modulation in direct DM detection 



has recently received further impetus from the CoGeNT 
experiment, located in the Soudan Mine, MN, USA. The 
collaboration has published its analysis on a 458 day run 
with a 440 gram Ge detector [TT]. After removal of cos- 
mogenic background contamination, the data exhibits a 
seasonal variation peaking around the middle of April. 
As we shall see below the modulation is manifest only 
in the higher energy part of the recoil spectrum where 
the energy spectrum is rather flat with none of the usual 
features expected from the direct detection signal of DM, 
elastic or otherwise. The unexplained exponential rise in 
the recoil spectrum at lower energies, which has stirred 
quite a commotion in the recent past, shows no evidence 
of annual modulations. Making use of the time-stamped 
data provided by the collaboration we address the mod- 
ulation part of the signal (see also ref. ^22, ,23J for prior 
analyses) and investigate the potential role of cosmic ray 
muons for the CoGeNT detector. 

In light of the significant advances in sensitivity from 
direct detection experiments achieved in the past decade 
and future improvements expected in the next one, we 
also address the question of whether the time-series of 
a signal may encode additional evidence that it is DM 
which is scattering on a target. In particular we point 
out that higher harmonics may be present in an annually 
modulated signal. We show how this signature manifests 
itself in the scattering rate and explore to some extent 
its sensitivity on the assumed halo parameters. 

The paper is organized as follows: We start with sec- 
tion [n] by establishing statistical evidence for periodic 
variation in the data sets under consideration. Restrict- 
ing ourselves to a sinusoidal form of the signals, in sec- 
tion [TTl] we allow both phase and period to vary and ex- 
plore the inferred values of these two parameters from the 
different datasets. In section |IV| we drop the assump- 
tion of a sinusoidal form and instead directly compare 
the time series of the various data sets by performing 
a correlation analysis. Section |V] discusses the effect of 
higher harmonic modulations as a new diagnostic tool 
when searching for dark matter in the lab. We summa- 
rize the main conclusions and findings of the paper in 
section IVII 



II. REJECTING THE NULL HYPOTHESIS 

The existence of a periodic signal in the DAMA and 
LVD data sets is clear even without any sophisticated sta- 
tistical analysis. The situation is somewhat more subtle 
in the case of CoGeNT. Both for the purpose of com- 
pleteness, and because it reveals interesting differences 
between the data sets, we begin our analysis by testing 
the null hypothesis of no signal (i.e. only pure noise) in 
each of the data sets. The most convenient way of doing 
such a spectral analysis is by looking at a periodogram 
of the data. 

The classical periodogram of Schuster [21] allows for 
a calculation of the power spectrum of a signal even in 
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the case of uneven sampling of the data. In ref. 
Scargle showed that a modified version of the classi- 
cal periodogram, which we call the Lomb-Scargle (LS) 
power spectrum (and sometimes just power spectrum), 
results in the same well defined statistical behavior as 
the Fourier power spectrum used in the case of even 
sampling. The LS power spectrum allows one to re- 
ject noise at a 1 — a confidence level for a single, pre- 
selected frequency by demanding a power level in ex- 
cess of Za = ln(a~^). If one is testing N independent 
frequencies then the power level required increases to 
Za = — In [l — (1 — a)-'^/'^] which contains the statisti- 
cal penalty for inspecting more than one frequency. This 
means that when testing for the null hypothesis over a 
range of frequencies we must employ some care in choos- 
ing the frequencies to be tested if we want a meaningful 
statistical interpretation. If the time series is evenly dis- 
tributed then the standard choice for the frequencies in 
the discrete Fourier transform results in independent fre- 
quencies. Since the data sets we consider are not grossly 
uneven and approximately cover an integer number of 
years we chose w„ = tuof with the fundamental frequency 
ijjp = 27r/T where T is roughly the range of years covered 
by the data sel[^ We checked that this results in low cor- 
relation among the test frequencies using the procedure 
described in App. D of Ref. pS]. 
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FIG. 2. The Lomb-Scargle power spectrum of the LVD data 
as a function of the period (solid- blue). Also drawn are the 
99% confidence lines for excluding the noise hypothesis for a 
single frequency (lower) or any of the frequencies shown (top) . 
Note the substantial power present for modes with a period 
greater than a year. The dashed-red curve shows the efiect of 
subtracting the yearly mean from the data on a yearly basis, 
as is done by the DAMA collaboration. The subtraction has 
little effect on the high frequency modes, but results in a 
strong damping of the long period modes as can be expected. 



A. DAMA and LVD data 

We start by considering the measurements of the niuon 
flux by the LVD experiment |19\ The average inte- 
gral muon intensity underground is reported as (/^) ~ 
3 X 10~^m~^ sec~^ with an annual variation of ~ 2% in 
amplitude as can be seen from Fig. [T] The data spans 
a total of eight years with more than 2800 data points 
which we obtain by digitizing Fig. 2 of [T^. 

The solid line in Fig. [2] shows the power spectrum ob- 
tained from the full LVD data set. Aside from the yearly 
modulation, which is plainly visible from the time se- 
ries itself, the power spectrum also speaks unequivocally 
of the existence of temporal variations on time-scales 
greater than a year. That same figure (dashed line) shows 
the effect of subtracting the mean intensity from the data 
on a yearly basis as done by the DAMA collaboration 
with their own data. The figure makes it clear that such 
a procedure would mask most of the power at periods 
much greater than a year. Nevertheless, we note that 
substantial power remains in modes with periods a little 
over one year. 



Figure |3] shows the power spectrum for the 
DAM A/LIBRA 2-4 keVee data set. In contrast to the 
LVD spectrum, there is little power at periods greater 
than a year. This, however, may simply be an artifact of 
the way the DAMA/LIBRA modulation data is obtained. 
The DAMA collaboration calculates the residuals by sub- 
tracting the mean on a yearly basis for each cycle. As 
mentioned in the previous paragraph such a procedure 
tends to dampen power at periods much greater than a 
year. Nevertheless, the absence of any power even for pe- 
riods slightly over a year, as seen in the LVD data above, 
is interesting and serves as the first distinguishing feature 
between the two data sets. Unfortunately, it is difficult to 
assign a quantitative significance to this difference with- 
out the full, unsubtracted time series. This issue further 
motivates the release of the unsubtracted data by the 
DAMA collaboration to allow for a proper comparison of 
the power spectrum. Finally, we note that no biannual 
mode (T = 1/2 year) is present in the DAMA/LIBRA 
power spectrum. We will comment on this further in 
Sec. El below. 



B. CoGeNT data 



^ The power spectra in this section all display a peak at a period 
of exactly one year. We caution the reader that this does not 
imply the best fit value would be exactly one year, but is simply 
an artifact of the frequencies we chose to sample when looking 
to reject the null hypothesis. 



The CoGeNT collaboration has released the time- 
stamped data of their 442 live-day run [TT] for public 
use. This makes the computation of the LS diagram in 
principle straight-forward. However, the data also suf- 
fers from background activity of electron capture decays 
of cosmogenically activated long-lived isotopes. A mea- 
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FIG. 3. The Lomb-Scargle power spectrum of the DAMA- 
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CoGeNT data as a function of the period for the three energy 
regions: 1.6 — 3.0 keVoc (solid-black), 0.9—1.6 keVco (dashed- 
blue), and 0.5 — 0.9 keVcc (dotted-purple). 



surement of higher energetic K-sheh captures together 
with reported ratios of L- to K-shell decays allows one 
to correct for the L-shell activity in the DM acceptance 
region. 

Based on the expected cosmogenic activity in the Co- 
GeNT dataset it seems reasonable to divide the low en- 
ergy data into three regions: 0.5 — 0.9 keVee, 0.9 — 
1.6 keVeo, and 1.6 — 3.0 keVee- Considerable cosmogenic 
activity is observed in the middle region. In contrast only 
a small amount of cosmogenic activity is expected in the 
low region (0.5 — 0.9 keVco) and very little if any is ex- 
pected in the high region (1.6 — 3.0 keVce)- Since the par- 
ticipating isotopes are rather long-lived {ti/2 ^ 200 days) 
they are expected to result in substantial power at long 
periods in the corresponding Fourier power spectrum. 
We have verified that this is indeed the case. However, 
since we are interested in more localized phenomena in 
frequency space, such as annual modulations, we need to 
subtract the cosmogenic activity from the data. We have 
done this on a daily basis based on the reported activ- 
ity levels We have verified that the results remain 
qualitatively the same even when adopting a different 
subtraction procedure, where we fit for the exponentially 
decaying component. 

Figure |4] shows the power spectrum for the subtracted 
data in the three energy regions. The power spectrum 
was calculated in 30 equally spaced frequencies with a 
fundamental frequency of 2 years. The null hypothesis 
of noise can only be confidently rejected for the 1.6 — 
3.0 keVeo region, where considerable power is present at 
a period of about one year. However, the same cannot be 
said about the lower energy regions where no significant 
power is observed. Thus, one can claim the detection of 
annual modulations at about 99% confidence level only 
in the 1.6 — 3.0 keVcc region. 



III. PHASE ANALYSIS 

Aside from the period, the second most important 
characteristic of the oscillations observed by the DAMA 
experiment is the phase of the signal. Dark matter colli- 
sion rate with the detector is expected to peak on June 
2"*^, corresponding to to = 152 days after January I''*. 
In this section we investigate the phase associated with 
the oscillations seen in the DAMA and CoGeNT data 
and compare them to the phase of the muon intensity 
modulations. This comparison was already attempted 
by the DAMA collaboration itself, however, it was criti- 
cized by refs. [171 HI] on two accounts. First, in their fit 
to the data the DAMA collaboration fix the period and 
allow only the phase to float. Second, the underlying 
signal may not be truly sinusoidal which may invalidate 
the statistical inference of a fit to a sine function. We 
address the first issue in this section and investigate the 
second problem in the next. 



A. Prequentist Approach 

We start with a simple frequentist approach and inves- 
tigate the effect of departing from strict annual period- 
icity (i.e T is not necessarily 365 days) on the DAMA 
phase to. Under the premise of a sinusoidal signal, 
A X cos[uj{t — to)] with uj = 2tt/T, we fit amplitude A 
to the DAMA/LIBRA residuals by minimizing the usual 
function while scanning over period T and phase to ■ A 
confidence region in to and T can be constructed from the 
profile likelihood method |27| which effectively removes 
the dependence on the nuisance parameter A. This way, 
we first obtain the global minimum Xmin from which the 
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FIG. 5. Confidence regions in period T and phase to ob- 
tained from a fit to tfie DAMA/LIBRA residuals. Tlie 
various lines as labeled correspond to the three different en- 
ergy binnings provided by the DAMA collaboration with the 
time origin set to Jan 1, 2003. The shaded region illustrates 
the effect of shifting the time origin to Jan 1, 2011 for the 
2-4 keV residuals and the light dot-dashed ellipse shows the 
time origin shift to Jan 1, 1995. 

confidence region is obtained by requiring 

x2(i,io,T) <xL„ + Ax^ (1) 

where A is the maximum hkeUhood estimate for A at 
each point (to, T^- For a coverage probabihty of 95% one 
chooses Ax^ = 5.99. Figure [s] shows the result of fitting 
the DAMA/LIBRA residuals in the (T, to)-plane. The 
three ellipses give the 95% C.L. regions for the various 
energy binnings as provided by the DAMA collaboration. 
As can be seen the data — at the required confidence — are 
not necessarily consistent with the dark matter interpre- 
tation (to = 152 days and T = 1 yr as indicated by the 
thin gray lines). As we argue in the next paragraph, the 
interpretation of such fits have to be handled with some 
care. 

When allowing the period to float it is important to 
realize that statements about the inferred phase become 
dependent on the starting date. First, there is the ob- 
vious effect that the phase is measured in days with re- 
spect to Jan 1^', but if the period is not one full year, 
then which year is used as the origin affects the phase. 
Comparison of the phase in days between experiments 
with different origins is meaningless until the origins are 
made to coincide. This effect is simple to correct for and 
requires that we use the same time origin for the different 
data sets. In Fig. [5] the chosen time-origin for the three 
ellipses showing the various energy bins is Jan 1, 2003. 



There is another, more subtle issue that affects the 
determination of the phase when the period is allowed 
to float. As can be seen from the ellipses in Fig. [5] the 
DAMA data exhibits an anti-correlation between the pe- 
riod and the phase. Fits with periods longer than a year 
strongly favour phases smaller than to ^ 150 days and 
vice-versa. Note, however, that the sign of the correla- 
tion itself depends on the chosen origin. When shifting 
the latter forward such that the DAMA/LIBRA data lies 
in the past, to and T become positively correlated in- 
stead. This is illustrated in Fig. [5] by the filled ellipse 
obtained from the 2-4 keV energy bins for which we have 
chosen Jan 1, 2011 as the origin. It is important to note 
that the distribution in to, given a certain period, de- 
pends on the chosen time origin. Ideally, we would like to 
make coordinate-independent statements and care must 
be taken when interpreting the results. This issue further 
motivates the correlation study provided in the next sec- 
tion. 

Given the above, before proceeding and comparing the 
LVD and DAMA data sets we must decide on a common 
time origin. The starting date used by the DAMA collab- 
oration (Jan I''*, 1995) is six years apart from the LVD 
data (Jan 1**', 2001). In our analysis we will concentrate 
on the DAMA/LIBRA data, which started on September 
9"^, 2003, and so a sensible choice is January 1*^*, 2003 as 
the origin since it has sufficient overlap with both data 
sets. To quantify the level of agreement of the muon- 
induced background hypothesis with DAMA we switch 
now to a Bayesian approach. This allows for a convenient 
generalization of the Lomb-Scargle periodogram into the 
two dimensional phase-frequency space as discussed in 
AppendixjA]. Such an approach is particularly convenient 
when one wishes to incorporate further assumptions on 
the provided data sets (e.g. choosing priors on the pe- 
riod, phase, and amplitude). We will not incorporate any 
such assumptions in the current analysis, since we would 
like to keep it as unbiased as possible. Moreover, we have 
verified that the conclusions arrived at below remain the 
same when we employ a frequentist analysis instead. 



B. Bayesian Approach 

The Bayesian approach can conveniently generalize the 
Lomb-Scargle periodogram and allow for the testing of 
different hypotheses and priors. In this section we re- 
strict ourselves to modelling the data with a simple si- 
nusoidal function /(t) = Ay. cos [ti;(t — to)] but following 
through the steps presented in Appendix |Aj a general- 
ization to more complicated test-functions is in principle 
straightforward. 

We are mainly interested in the relationship between 
phase to and period T. The posterior probability 
P({a;, to}||d), i.e. the probability of observing frequency 
a; = 2-K jT ^ and phase to given the data d is obtained 
by integrating the full posterior P({A, w, to}| |(i) over the 
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amplitude A restricted to positive value^ 
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p = ^cos^uj {ti - to) , 

■i 

h = di cos uj {ti — to) . 



(2) 

(3) 
(4) 



Before proceeding and using ([2| in the fuU two- 
dimensional period-phase space we can make contact 
with the procedure which is usually employed, i.e. fix- 
ing the period to one year, T = 365 days. Thereby we 
impose a delta function prior on the period centered at 
one year from which we consequently obtain the poste- 
rior probability in phase only, P({io}|M)- We find that 
the DAM A/LIBRA and LVD data do not agree with re- 
spective values 



to(DAMA) = (13f ± 13) days, 
to (LVD) = (187 ±2) days. 



(5) 
(6) 



The tight error bars on the LVD data may at first seem 
surprising. To better motivate this number we consider 
that given the 8 years of data, one can easily determine 
by eye a phase shift of about 20 days. Statistics allows 
a further reduction of about an order of magnitude be- 
cause of the very large number of points available in the 
dataset. However, it should be kept in mind that the un- 
certainty quoted is only statistical and does not take into 
account the possible systematic issues associated with the 
obvious presence of temporal variations on larger time 
scale^ In fact, using the error bars of the LVD data, 
the resulting of a fit to a sinusoidal function is ex- 
tremely poor. Hence, the numbers quoted above should 
really only be understood as the best inference on the 
underlying parameters of a model consisting of a single 
frequency. 

To see if the discrepancy between DAM A/LIBRA and 
LVD remains when relaxing the assumption of a strict 
annual periodicity, we now evaluate ([2| for the two data 
sets. The result is shown in Fig. [6] This figure clearly 
shows that the inferred modulations seen in the two data 



It is possible to obtain a more compact expression for the poste- 
rior without the error function by integrating over all amplitude 
values both positive and negative. However, that leads to a ±tt 
ambiguity in the phase. While this degeneracy is easily remov- 
able by inspection, we prefer to avoid this complication in what 
follows. 

Moreover, there is a systematic error associated with the digiti- 
zation of the data that amounts to an additional uncertainty of 
a couple of days. 
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FIG. 6. A comparison of the period-phase posterior for 
DAMA/LIBRA (bottom) and the LVD phase (top). Shown 
are the 68%, 90%, and 99% credibility ellipses. 



sets — if interpreted as sinusoidal modulations — are in- 
compatible with each other. 

A similar analysis can be done for the CoGeNT datsQ 
as shown in Fig. [7j where we chose to use the entire low 
energy range (0.5 — 3.0 keVce). As discussed in the pre- 
vious section, evidence for modulation is only present in 
the restricted range of 1.6 — 3.0 keVoo. However, iso- 
lated power on an annual time scale is present in the 
entire range and it is therefore not unreasonable to em- 
ploy the full range when looking to make inferences about 
the modulation parameters. CoGeNT's credibility el- 
lipses are compared with those obtained from the MINOS 
data [29] . From there, we arrive at similar conclusions as 
before, namely, that the muon and direct detection data 
sets seem incompatible We have verified that when con- 
sidering the restricted range 1.6 — 3.0 keVco in the Co- 
GeNT data we arrive at the same conclusions. We note 
in passing that the orientation of the MINOS ellipses — 
which are obtained with a time origin of Jan I''*, 2010 — is 
similar to the one found in Fig. [5] for DAMA once one 
shifts the time origin to the future of the data set. 

The credibility contours in Fig. [7] indicate a range of 
viable parameters which is almost twice as large as what 
has been quoted in the CoGeNT release paper [TT] . For a 
direct comparison we thus also perform a frequentist fit to 
a cosine function (plus a constant.) For example, when 
using the Minuit package [30j for the ^^-minimization. 



* A Bayesian approach to the CoGeNT data has very recently also 
been taken by |28l . 
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FIG. 7. The period-phase posterior for the full low energy 
CoGeNT data (0.5 - 3.0 keVee) [H] and MINOS data [H]. 
Shown are the 68%, and 90% credibility ellipses. The time 
origin in both cases was chosen to be January 1^^, 2010. 



ertheless, we now set to investigate this issue in a way 
that does not rely on assuming any particular functional 
form for the underlying temporal variations. 

Moreover, even if the timing between muons and 
DAMA (and possibly CoGeNT) are incommensurate at 
first sight, could it be that the underlying stochas- 
tic nature of the background process alleviates the ob- 
served tension in the annual phase? This explanation for 
DAMA was suggested recently in an interesting paper by 
Blum 121) . Assuming a simple, generic model for how 
a muon-sourced background may be realized, the distri- 
bution in phase originating from the Poissonian process 
is claimed to be compatible with the one observed by 
DAMA. 

In this section we address the above issues by per- 
forming a correlation analysis. In particular we follow 
Ref. [21] and generate mock data for DAMA based on 
the muon hypothesis. We show that the resulting mock 
data exhibits a substantially higher degree of correlation 
with the actual muon data than does the real data from 
DAMA thus ruling out the hypothesis. We then go on 
and use a similar correlation analysis to show that the 
muon hypothesis is also an unlikely cause for the modu- 
lations seen by CoGeNT. 

A. DAMA 



we obtain T = 350 ± 26 days, = 110 ± 13 days, and 
A = 17% ± 4% for a time binning resulting in 21 bins. 
This is in good agreement with [TT]. However, it should 
be stressed that the quoted errors are obtained from a 
default value of Ax^ = 1 which does not reflect the 
increased freedom of fitting more than one parameter. 
Furthermore does not grow very large with respect to 
its rather small minimum value, x^^^/ d.o.f . = 6.8/17 
and even a fit to a constant rate yields a reasonable 
Xmin/^-O-/- = 26.9/20. This can be traced back to large 
error bars in conjunction with limited data. As a re- 
sult of this, coirtours with nominally larger confidence 
will rapidly open up the parameter space. This, how- 
ever, does not signal real compatibility of the data sets 
because of the aforementioned reasons. 



IV. CORRELATION ANALYSIS 

A valid concern with regard to the analysis above is 
the underlying assumption of sinusoidal variation. The 
power-spectrum of the LVD data. Fig. [2] makes it clear 
that power exists in modes with periods larger as well 
as smaller than one year. One may then worry that the 
phase comparison discussed in the previous section suf- 
fers from a systematic misinterpretation. We do not en- 
tirely endorse this concern because it is difficult to under- 
stand how the very prominent annual modulation in the 
LVD data, with its very definite phase, can be masked by 
the much weaker modulations at other frequencies. Nev- 



Since there is considerable overlap between the 
DAM A/LIBRA and LVD data, it is straightforward to 
define Pearson's coefficient of correlation. We first bin 
the LVD data according to the DAMA/LIBRA bins. 
Then the sample's correlation coefficient is defined as, 

1 Ellr {X, ~ X) (y. - Y) 
rx,Y = 7 [') 



1 



axo-y 



where n is the number of overlap bins, X and Y are 
the samples' mean, and ax,Y are the samples' standard 
deviations. The first statistical question we address is 
whether we can exclude the no-correlation hypothesis. 
The answer to this question lies in the statistic. 



t 



r\fn 



(8) 



which — in the case of the null-hypothesis — is known to 
follow the Student's distribution with d.o.f — n~2. The 
DAMA/LIBRA and LVD data share 39 temporal bins 
and so the 90% (99%) two-sided exclusion limit on the 
statistic is |t I > 1.687(2.715). For the correlation between 
the two data sets we find a value. 



''LVD, DAM A = 0.44 



^LVD,DAMA = 2.95, (9) 



which allows us to exclude the no-correlation hypothesis 
with confidence level greater than 99%. This, however, 
should come as no surprise since both datasets exhibit 
strong annual modulations with a similar phase. Any 
such samples will exhibit a correlation at some level. 
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FIG. 8. The population distribution of the phase for the 
model by Blum, Eq. ( |10[ ) compared with the best fit to data 
obtained by the DAMA collaboration, tg — 136 days, shown 
as the red vertical line. 



but that of course does not imply that they are indeed 
causally related. The more interesting question we would 
like to address now is whether we can exclude causation. 
This question can only be answered in a model-dependent 
way. 

The model we consider was presented by Blum in 
Ref. [H]. It assumes that the DAMA count rate Si in 
a time bin of width Ati in an energy range AE for a de- 
tector of mass Al is related to the muon flux /^.j during 
that time by. 



yN^ 



where N,, 



' MAEe.Ati ' 
is Poisson distributed with mean, 



A-cffIfj..ieiAti 



(10) 



(11) 



Here y is the number of signal counts/muon (yield), Acs 
is an effective area in which muons arc collected, and 
Ci ~ 60% is the duty cycle during time bin ti. As was 
argued in |21[ , direct hits of muons in the detector would 
require y w 50 and would lead to a spread cr, /si in events 
which is about a factor of five larger than what is actu- 
ally observed. Since this is in clear conflict with the data, 
the other possibility that we are going to consider is the 
secondary effect muons can have due to the spallation 
reactions on nuclei and the neutron flux they induce. In 
this case muons can be collected in a larger area, say, 
Acs « 10 m^, which requires an order one yield y ~ 2 
only. Using the latter numbers we generate 10^ realiza- 
tions of DAMA data based on this model. 

We first attempt to make contact with [2T] by plot- 
ting in Fig. [s] the distribution of to obtained from a x^- 
fit to a sinusoidal function with floating period, phase, 
and amplitude. There is little resemblance with the 
equivalent Fig. 3 presented in [2l]. The latter shows 
a very broad distribution with substantial support be- 
tween 90 < to/days < 250 and a peak at to ~ 100 days. 
We find from our Fig. [8] that to is normally distributed 
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FIG. 9. The population distribution of the Fisher Z-transform 



for the model by Blum, Eq. ( 10 1 as compared to the sample 



value of Z = 0.47. Based on the population distribution, the 
model can be excluded with a confidence exceeding 99%. 



with a mean of (to) — 187 days and a = 27 days. We 
believe that the disagreement with Ref. [21] is due to a 
different choice of the time-origin. As we discussed in 



section III once we allow the period to float, the phase 
to loses its absolute meaning and the allowed region in to 
becomes sensitive to the origin of time. By shifting the 
origin from Jan I''' 2003 to the one DAMA uses when 
quoting their data, Jan 1^* 1995, we find that the distri- 
bution in to indeed broadens significantly similar to the 
one found in |21| . Though it may certainly be better to 
choose the time origin in 2003, the discrepancy signals 
a more serious conflict: the distribution in to does not 
provide us with a robust test statistic when assessing the 
compatibility of the mock data with the observed phe- 
nomena. 

A better way to assess the (in)compatibility between 
the muon flux and DAMA is to use the generated set of 
realizations and evaluate the degree of correlation itself. 
This way we remain completely coordinate independent 
in our statements and we can attempt to exclude the 
model by showing that it implies a level of correlation 
much higher than what is actually implied by the data. 
This is precisely what we now labour to show. In general, 
if the model allows an exact calculation of the population 
correlation coefficient vq then the Fisher Z transform, 



ilog 



1 



1 - r 



(12) 



is approximately normally distributed with mean and 
standard deviation given by. 



llog 



1 + ro 
l-r-o 



(13) 



We will not use these results below since the model we 
consider does not easily allow for an analytic evaluation 
of the population correlation coefhcient tq. Instead we 
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will utilize the numerical realizations of the model and 
calculate the distribution explicitly. As it turns out, how- 
ever, the numerical results nicely agree with the theoret- 
ical expectations for the shape of the distribution. 

Figure [9] shows the results of the computation of the 
Fisher Z-transform for each of these realizations. The 
distribution is compared to the sample's Z-transform 
Z = 0.47 obtained from (|9]). The numerical realizations 
reveal a much larger expected correlation than what is 
observed in the data. Therefore the model is excluded 
with a confidence greater than 99%. 



B. CoGeNT 

We now move on to investigate the hypothesis that 
the reported CoGeNT signal is caused by cosmic muons. 
In this respect, it is important to note that the apparent 
modulation fraction of the CoGeNT signal is most promi- 
nent in the high-energy bin 1.6 — S.OkeVco with a value in 
excess of 10% [TI]|22]- Seasonal muon-flux modulations 
of that order have indeed been reported by IceCube [3T] . 
However, the (somewhat) milder climatic conditions at 
the Minnesota Soudan Mine location exhibit a variation 
of at most 4% |29j. Therefore, unless the scaling of the 
signal is — for some unknown reason — non-linear in inci- 
dent muon flux or unless the performance of the detector 
was not stable throughout the data taking period, it is 
not possible to generate, say, a 16% CoGeNT modulation 
from a lesser modulated sourcing process. This state- 
ment is independent of the potential presence of further 
background. Even though this argument makes a muon 
explanation of CoGeNT rather unlikely, we can still pro- 
ceed and look for a quantitative answer based on the 
temporal structure alone. 

The main complication associated with such an anal- 
ysis is the fact that although measurements of the un- 
derground muon flux are available from the MINOS ex- 
periment pQj, the data has no temporal overlap with 
CoGeNT. We circumvent this issue by relying on nearby 
atmospheric temperature data rather than on measure- 
ments of the muon flux itself. As is well known, the 
underground muon flux is tightly correlated with the 
(stratospheric) temperature. We therefore choose to 
directly evaluate Eq. ([t]) between the effective atmo- 
spheric temperature parameter Tes and the background 
subtracted CoGeNT data. We refer the reader to Ap- 
pendix [B] for further details on this procedure. 

A correlation analysis is particularly appropriate in 
this case since CoGeNT's relatively short data taking pe- 
riod (458 days) naturally limits the significance of any 
statements about the annual periodicity of the signal 
as evinced in Fig. [7j However, the correlation analysis 
we employ requires binning the data, a procedure that 
is complicated by the need to subtract the cosmogenic 
background which is responsible for the majority of the 
event rate in the 0.9 — 1.6 keVoc region of the CoGeNT 
data. One may worry that such subtractions introduce a 
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FIG. 10. Coefficent of correlation r between Tch and CoGeNT 
data as a function of bin size in live days. The various fluctu- 
ating lines correspond to different binnings in energy. None 
of the chosen energy regions show any significant degree of 
systematic, positive correlation. 



bias depenedent on the choice of binning or energy range. 
Moreover, the unmodulated component of the subtracted 
CoGeNT data exhibits a rise towards low energies that 
is not shared by the modulated part, which if anything 
seems larger at higher energies. Thus, it is not even clear 
which part of the energy spectrum could be a result of 
the muon flux. In light of these complications we con- 
sider various time- and energy binning in the hope of 
thoroughly covering the different possibilities. 



Figure [TO] shows the correlation coefficient as a func- 
tion of bin size in live days for various energy ranges. If 
any (positive) correlation was present, which would cor- 
roborate the hypothesis that muons are responsible for 
the observed signal, one would expect significant degree 
of correlation independent of the bin size. With the co- 
efficient of correlation fluctuating near zero, the opposite 
is observed. The step-like red lines delineate the (two- 
sided) 90% CL for rejecting the null-hypothesis. The 
data is therefore perfectly consistent with the null hy- 
pothesis of no correlation. For larger bin sizes r begins 
to fluctuate more strongly. Although it is easier to reach 
some degree of correlation with fewer bins it correspond- 
ingly becomes harder to reach a given level of significance 
which is indicated by the opening of the red lines. 
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V. BIANNUAL MODULATIONS 

As discussed throughout this work, one of the hall- 
marks of dark matter direct detection is the annual mod- 
ulation in the recoil rate, which is the result of the Earth's 
motion with and against the WIMP wind 13 [5] . How- 
ever, as we show in this section, one also expects higher 
harmonics (e.g. a period of half a year) to be present 
at some level and those may prove useful as additional 
confidence builders if a signal is detected. 



A. Harmonic Analysis 

The differential recoil rate, dR/dE^^ is a function of 
the Earth's velocity relative to the rest frame of the dark 
matter halo, which using the notation of Ref. [31| is ap- 
proximately given by, 



Vobs 



V0 + V® (ei cosw (t - ti) + 62 sinw {t ~ ti)) (14) 



Here Vq is the velocity of the sun relative to the halo, and 
ii (£2) is the velocity unit vector of the Earth's motion at 
ti = March 21"' {ti -I- year/4). Due to Earth's orbit, the 
second term proportional to the Earth's relative velocity 
to the sun, ~ 30 km/s, modulates with frequency 
If the dark matter velocity distribution is 



271- 
1 year 



isotropic, then the dark matter speed distribution in the 
earth's frame, f{v), depends solely upon the magnitude 
of Vobs- To a good approximation the magnitude of the 
velocity is given by 

|vobs| « IvqI + ^V©cosw(t-to). (15) 

Thus, dark matter scattering rate is a periodic function 
with a fundamental period of 1 year, but like any periodic 
function it may contain higher harmonic modulations of 
any frequency a;„ — nuj. We therefore expect modula- 
tions in the scattering rate with periods of 1/2, 1/3, ... a 
year corresponding to biannual, triannual and higher fre- 
quency modulations. In addition, there are harmonic cor- 
rections to I Vobs I itself. They arise predominantly from 



corrections to ( 14 1 due to smaller effects such as the ellip- 
ticity of the earth's motion around the sun. Once these 
are taken into account, they result in observable phase 
shifts between the different harmonics. 

We can learn more about the higher harmonics, by 
expanding in the parameter e^, = V^/2v(7) ~ 0.06. The 
time dependence of scattering arises solely due to the 
velocity dependence through 



dR 

oc 



dE 



f{v) 



dv 



cosa;(t - to)]" (16) 

=o,i,--- 

^ c„(-i;,nin) cos [na;(i - to)] 



Here we have used the fact that trigonometric identities 
allow powers of cos [uj{t — to)] to be re-expressed as sums 
of cos [nw(t — to)] - In EQ- ( |16| ) ^'min is the minimum dark 
matter speed in the lab which can deposit a recoil energy 



V - 1 ( ^nEr ^ ^ 



(17) 



for a nuclear target of mass mjv and dark matter-nuclei 



reduced mass fj^Nx = rriNm^/ {mi^ + m^)\ in (17 1 we 



have included the possibility of inelastic scattering with 
a splitting i5 in the dark matter sector [5J [S] . The com- 
putation of the harmonic coefficients requires a choice of 
velocity distribution, and in what follows we consider the 
distribution proposed in [M] . 



fk{v) oc 



exp 



kvl 



- 1 



Q{v,,,-v) (18) 



Here, Q{x) is the Heaviside step function, and we con- 
sider the dispersion i;o = 220 km/s [55H57] and the es- 
cape velocity Vesc = 600 km/s [55] . 



In Fig. 11 (top), we illustrate the behavior of the har- 
monic coefficients, c„, for the fc=1.5 velocity distribution. 
In order to plot the amplitudes on a logarithmic scale, we 
have taken the absolute values of c„. The troughs in the 
plots indicate when c„ is changing sign. At high enough 
minimum velocity, Umin, the harmonic coefficients are all 
positive. Interestingly, there are n sign changes for c„, 
which can be understood from the behavior of the veloc- 
ity distribution f{v). 



In the approximation of Eq. ( 16 ), the phases of all the 
higher harmonics are the same as the annual modula- 
tion. However, as mentioned above, there are higher or- 
der corrections to the magnitude of the observer's veloc- 
ity, Eq. ( 15 ), that allow the phases of the higher harmon- 



n=0,l, 



ics to deviate from the phase of the annual modulation. 
Indeed, using an accurate parameterization of the earth's 
velocity |33j . we find appreciable temporal shifts of the 
biannual and triannual mode with respect to the annual 
one. This is shown in the middle panel of Fig. [TT] The 
abrupt changes seen in the figure indicate sign changes 
in the harmonic coefficients c„, which can be thought of 
as phase shifts by 1 yr/2n. The ambiguities in the phase 
shifts are fixed by choosing the values which are closest to 
the annual one. These phase shifts serve as an additional 
signature which can be used to discriminate a potentially 
positive signal from other sources. 

In the bottom plot of Fig. [TT] we plot the modula- 
tion amplitude ratio c„/ci for the biannual and triannual 
modulation divided by the annual modulation. There are 
two regions where the higher harmonics are as important 
as the annual modulation. The first region occurs near 
200 km/s around the zero in the annual modulation co- 
efficient ci. Here, the biannual amplitude C2/C0 is at 
the per mille level which is small but potentially observ- 
able with a large amount of exposure. The second region 
where higher harmonics are important is at high veloci- 
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ties, above the escape velocity cutoff. At these high ve- 
locities, there are substantially less dark matter particles 
in the winter than in the summer. Thus, the scatter- 
ing rate can vanish entirely in winter, with a chopped 
cosine behavior, which requires more harmonics to fully 
describe the approach to zero. Although this would be 
obvious in a experiment with no background, the pres- 
ence of background can mask these higher harmonics by 
eye, so it is useful to look for them. 

We have also looked at the behavior of these modu- 
lation amplitudes for the Standard Halo Model and also 
for a k = 3.5 distribution [33] . The behavior at low 
Vmin is vcry similar, since that region depends mostly on 
the velocity dispersion and not on the behavior near the 
tail. The high velocity behavior enhances (suppresses) 
the higher harmonics for sharper (shallower) escape ve- 
locity cutoffs as seen in the k = 3.5 distribution (Stan- 
dard Halo Model). 

A detailed analysis of the detectability of these higher 
harmonics and their ability to be faked by backgrounds 
is beyond the scope of this paper. Here, we restrict our- 
selves to a few remarks. Any background that modu- 
lates annually, with a small modulation amplitude will, 
through the same Taylor expansion argument, also have 
higher harmonics (although not necessarily all with sim- 
ilar phase). What is peculiar about dark matter is that 
the higher harmonic amplitudes are enhanced at high 
'^inin due to the escape velocity physics. This could be 
mimicked by background only if there was some effect 
during winter which severely suppressed the background 
appearing as a nuclear scattering event. 

On the theory side, models with inelastic dark matter 
scattering lead to enhanced higher harmonic amplitudes 
as compared with standard elastic dark matter. A split- 
ting 6 ^ 100 keV leads to values of Umin which can be 
near the cutoff (see Eq. 17 1. In fact, we find that the 



biannual amplitude can be as large as 30% of the annual 
modulation amplitude for the 2-4 keVee bins of DAMA. 
Looking at Fig. [3| we see that this is consistent with the 
amount of biannual power seen in these bins, since the 
power ratio between biannual and annual modulation is 
(c2/ci)^ ^ 0.1. With a large increase in exposure, by a 
factor of 0(10 - 100), DAMA would have enough statis- 
tics to begin testing this and perhaps even see evidence 
for a nonzero biannual amplitude as we show quantita- 
tively in the next subsection. 

As a final remark, we note that it is possible to pre- 
dict the energy bin where the annual modulation should 
be suppressed, due to its change of sign, which would 
be an interesting bin to look for biannual modulation. 
For heavy dark matter that scatters elastically, the de- 
pendence on the dark matter mass drops out of the re- 
duced mass, so we can use Eq. [iTj to show that the 
energy bin Er = 90 keV(mAr/100 GeV) is where the 
annual modulation should be suppressed. Such high en- 
ergies might be hard to get substantial dark matter rates 
though, since the nuclear form factors tend to suppress 
such high energy recoils. The situation gets better for 
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FIG. 11. The velocity dependence of dark matter scatter- 
ing as a function of «niin for the A;=1.5 velocity distribution 
proposed in 04], with dispersion «o = 220 km/s and escape 
velocity Wesc = 600 km/s. The top plot shows the |cto| val- 
ues for n=0, 1,2,3 (top to bottom) where c„ is in units of 
(km/s)~^, the middle plot is the phase shift of the higher 
harmonics relative to the annual modulation, while the bot- 
tom plot shows the modulation amplitudes |c„/ci| for n=2,3 
(top to bottom). The horizontal lines in the bottom panel 
show the 90% C.L. upper limits on a hi- and triannual signal 
in the DAMA data. 



dark matter much lighter than the nucleus, as the energy 
bin is Er = 0.9 keV(TOx/10 GeV) ^(to at/ 100 GeV)-i. 
This requires knowledge of the dark matter mass to pre- 
dict, but experimentally it might be verified that the sign 
of the annual modulation amplitude changes sign in a cer- 
tain bin, thus it would be an interesting foUowup to see if 
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there is a biannual mode in that bin. In fact, the DAM A 
annual modulation amplitude is smaller in the first bin 
2 — 2.5 keVoc bin, so it would be interesting to see if there 
is any biannual modulation there. 



B. Testing the Signal Hypothesis 

In the presence of a signal, the probability distribution 
function of the power P in a particular frequency of the 
Lomb-Scargle power spectrum changes to 25, 39 , 

p(P;P,) = /o(2/PP:)cxp(-P-F,). (19) 

This formula assumes a normalized signal with power 
Ps that is measured in units of the variance of the noise. 
Given an observation of power Pobs we can reject a signal 
at a level of 1 — a by requiring that, 



Pr(P > Fobs) = 1 - 



p{P;Ps)dP = l-a (20) 



So, for example, from Figure [3] the DAMA/LIBRA data 
in the 2—4 keVee energy range has a power of Fobs — 0.57 
(1.8) in the biannual (triannual) mode. This implies a 
frequentist 90% CL upper limit of pbi^nnuai < 1.9(4.4) 
on the power of the biannual (triannual) mode of a puta- 
tive signal. Since DAMA claims detection of an annually 
modulating signal we can find the ratio of the higher 
harmonic power to the annual power. The power ob- 
served in the annual mode of the same energy bin is 11.0. 
Therefore the fractional power in the biannual or trian- 
nual mode must be lower than. 



^biannual ^ q -j^y^annual 
^triannual ^ q ^^annual 



at 90% CL. 



(21) 



Since the power is directly related to the square of the 
harmonic coefficient, this limit implies a corresponding 
limit on the harmonic coefficients, C2/C1 < 0.42 and 
C3/C1 < 0.63. These limits are shown in the bottom 
panel of Fig. |11| A similar analysis for the muon flux 
gives a 90% CL upper limit on the biannual mode of 
pbiannuai ^ g 3^ Figure^ Sincc the power observed 
in the annual cycle is PJ>-"""^ = 312 this implies a limit 
on the harmonic coefficients of C2/C1 < 0.17. 



VI. CONCLUSIONS 

This paper was devoted to a detailed time-series analy- 
sis of Dark Matter direct detection data as well as related 
datasets. The main findings are as follows: there is no 
evidence to support the idea that atmospheric muons are 
responsible for the signal that the DAMA collaboration 
observes; the annual modulations observed in the Co- 
GeNT data are statistically significant only in the higher 
part of the low energy spectrum (1.6 — 3.0 keVoo); no 
correlation of the latter with the expected muon fiux is 



observed; biannual modulations may serve as an addi- 
tional handle on any putative dark matter signal. 

With regards to the muon hypothesis discussed in 
refs. [ini HZl [3T] whereby atmospheric muons reaching 
the DAMA detector are responsible for the signal ob- 
served, our analysis indicates three significant difficulties 
with this idea: 

1. The power spectrum of the LVD data on the 
muon flux indicates considerably more power at 
periods longer than one year compared with the 
DAMA/LIBRA data. This remains true even after 
we follow the undesirable procedure of the DAMA 
collaboration whereby the residuals are computed 
by subtracting the mean of the data on a yearly 
basis. This procedure is undesirable since it masks 
power at long periods. We therefore urge the 
DAMA collaboration to release the unsubtracted 
data to allow for a full comparison of the power 
spectrum. Similar conclusions have been reached 
in |21]- 

2. When fitting the two datasets to a sinusoidal func- 
tion with variable amplitude, phase, and period the 
most likely regions for the two fits do not overlap. 
On this point we are in disagreement with ref. [2T\ 
which finds that the phase and period of the muon 
dataset can be in agreement with the DAMA sig- 
nal. We believe that the source of this discrepancy 
is the choice of time origin which becomes impor- 
tant once the period is allowed to float. 

3. A correlation analysis which is independent of any 
choice of functional form or flt to the datasets ex- 
cludes the simple, yet very general model presented 
by ref. [21j for connecting the muon flux with the 
DAMA signal. It indicates that such a model would 
imply a much stronger degree of correlation than 
what is in fact observed in the data. Thereby, 
this analysis also puts to question the hypothe- 
sis [121 [T7] itself that the variation of the muon 
flux has a casual connection with the one seen by 
DAMA. 

Aside from these three difficulties, which depend only 
on the temporal structure of the two signals, there also 
seems to be some difficulty reconciling the amplitude of 
the muon modulations with that of the DAMA signal. 
On the face of it, the amplitude of the modulations ob- 
served by LVD seems to agree well with the ~ 2% modu- 
lations as originally quoted by the DAMA collaboration. 
However, the full single-hit energy spectrum as measured 
in DAMA-LIBRA [5], seems to suffer from '"^K contam- 
ination in the signal region. Moreover, it was recently 
pointed out in ref. [ID] that the radioactive background 
rates in the DAMA experiment are likely larger than 
what was estimated by DAMA. If true, that would imply 
a smaller unmodulated amplitude of the putative signal 
and hence a larger modulation fraction. Thus, any model 
attempting to explain the DAMA signal as a result of the 
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muon flux will require some sort of nonlinearity in the 
signal yield. 

Given all of the above, we find the niuon hypothesis in 
its current form rather unpersuasive. Any further inves- 
tigations into this possibility should address the afore- 
mentioned difficulties. 

With regards to the CoGeNT data our analysis concen- 
trated on three different energy regions, 0.5 — 0.9 keVco, 
0.9 - 1.6 keVoc, and 1.6 - 3.0 keVeo. We find that the 
presence of annual modulations can only be firmly es- 
tablished in the 1.6 — 3.0 keVee region. The other two 
regions require a somewhat more careful treatment of the 
data to remove the contribution from the cosmogenic ac- 
tivity which contributes to the power spectrum at long 
periods due to the long-lived isotopes involved. Nev- 
ertheless, none of the procedures we have followed for 
subtracting the cosmogenic contributions resulted in a 
statistically significant yearly modulation signal in the 
subtracted data. The phase analysis of the modulations 
shows them to peak around the middle of April, which is 
somewhat earlier than expected from DM recoils, but the 
uncertainty on this number is still rather large. Indeed, 
we find that the confidence intervals in the full two di- 
mensional phase-period space are much too large for any 
significant statement to be made about the modulations 
observed. 

As in the case of DAMA we find no significant cor- 
relation between the CoGeNT data set and the yearly 
modulation of the muon flux. In the absence of published 
measurements of the latter for the time-period of the Co- 
GeNT run we used the atmospheric effective temperature 
inferred from climate data instead. The one-to-one rela- 
tion between the variation in temperature and muon-flux 
allows us again to indirectly disfavor the muon hypothe- 
sis as it applies to the modulations seen by CoGeNT. 

One can also wonder if muon induced backgrounds in- 
stead have a delayed signal, would that allow them to 
explain these dark matter anomalies? As shown in Ap- 
pendix [Cj production of long-lived radioactive isotopes 
does not result in any improvement. The largest allowed 
temporal shift of the signal peak is to the beginning of 
October. To complicate matters, the modulation ampli- 
tude of such a signal is always smaller than the muon 
amplitude. These limitations make it difficult to explain 
DAMA or CoGeNT with such muon backgrounds. 

Finally, we discussed a new diagnostic for dark matter 
direct detection experiments. Since the expected recoil 
rate of dark matter against matter in the lab is a general 
periodic function with a fundamental period of one year, 
it also contains biannual modulations as well as higher 
harmonics. The relative amplitude of the velocity mod- 
ulation is roughly V^/2vq ~ 0.06, which is rather small. 
The higher harmonics in the recoil rate have similar phase 
to the annual modulation, but are generally suppressed. 
However, we presented a few cases where they might be 
observed given sufficient exposure. At the very least, the 
biannual mode can be looked for and limits on its power 
can be obtained in a straightforward fashion. 
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Appendix A: Bayesian Version of the Lomb-Scargle 
Periodogram 

In this appendix we produce a derivation of the Lomb- 
Scargle periodogram based on Bayesian analysis. This 
result was first derived in Ref. [H]. We present it here 
in order to introduce the terminology, methodology, and 
notation used in deriving some of the other results, as 
well as in order for this paper to be as self-contained as 
possible. 

We begin by considering the oscillatory model func- 
tion / with a single harmonic frequency w, 



f{t) = A cos (ujt) + B sin (ujt) , 



(A-1) 



as the hypothesis for the signal. The data points are 
therefore given by. 



d^ - f{t,) 



(A-2) 



where the data is taken to have zero mean (d) = 0. The 
noise is denoted by rii which is assumed to be normally 
distributed with zero mean and variance a^. Under the 



hypothesis (A-1) the probability P{d\\{A,B,ijj}) of ob- 
serving the data d given the parameters {A, B, a;} is then 
directly proportional to the likelihood function 



N 



L{{A,B,u),d) ^\{^==c^Y> 



=1 ^ 



- f{u)f 



2a2 



1 



N 



cxp 



, \/27nT 

where N is the number of data points and 



_NQ\ 
2^2 ) 



(A-3) 



2 

N 



A di cos {ujti) + B di sin {toti 



1 

TV 



A^ ^ cos^ {ujti) + B^ ^ sin^ (wt^) 



2AB cos {ujti) sin {ujU) 



(A-4) 



Bayes' theorem then allows us to obtain the posterior 
distribution P ({w, A, i3}||d), i.e. the probability of the 
parameters {w. A, B} given the data d, 

P{{u:,A,B}\\d) cx P{{Lo,AM)L{{A,B,u:},d). (A-5) 

Assuming a fiat prior on the nuisance parameters A 
and B, we can integrate over A and B to obtain the pos- 
terior distribution of uj alone. This is a gaussian integral. 
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and would be trivial except for the cross term between 
A and B in Q. This term can be eliminated with the 
introduction of a phase-shift, r, by defining ti = — t 
and using the fact that {d) = 0, 



id' 
1 

N 



2 

N 



A di cos (wti) + B di sin (ojii) 

i i 

cos^ {uii) +B^Y^ sin^ (wt,) . (A-6) 



Effecting the gaussian integral over A and B gives, 
P{{uj}\\di) cx j dAdB P{A,B) L{{A,B,uj},d) 



1 



V2 



N-2 



1 



^^^^ ^E.cos2(^t,)E.sin2(u;t,) 
1 



2a' 

where the LS function is, 
LS(w) = ^ 



J2t cos2 [ujU 

1 



[iV(d2) - LS(a;)] 

cos (wti) 
sin (wti) 



(A-7) 



H 2 



(A- 



Y,, sin^ (wii 

The amplitude of the exponential in the posterior distri- 



bution, Eq. (A-7), is a rather mild function of the fre- 
quency and so the most likely value for the frequency is 
the one which maximizes the LS function. This result is 
the Bayesian version of the usual frequentist LS statistics. 

This approach lends itself to several generalizations 
and extensions. Ref. [41, derives the posterior for the fre- 
quency associated with a sinusoidal signal with general 
coefficients (not necessarily constants) . It is also possible 
to derive the posterior on the phase of a sinusoidal signal 
following similar steps to the above. 



Appendix B: Cosmic muon flux 

In this section we discuss the connection between cos- 
mic ray induced muon flux and its seasonal variation. 
Both, DAMA and CoGeNT, are operated deeply under- 
ground (> 2000 mwe) so that these detectors are well 
shielded from cosmic rays. Only muons with initial en- 
ergies ©(TeV) are able to reach the respective LNGS 
and Soudan Mine underground facilities. The highly 
energetic muons are sourced from primary meson (tt's, 
if's,. . . ) decays which themselves are produced when 
cosmic rays hit the atmosphere. Since the chance for pri- 
mary mesons to undergo further interactions grows with 
air density, there is an increasing competition between 
further degradation and decay when going to lower at- 
mospheric levels. 
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FIG. 12. Evolution of effective temperature parameter Tctt at 
the Soudan site as a function of time. The gray shaded regions 
indicate the overlaps with data-taking periods of MINOS and 
CoGeNT as labeled. 



The above chain of arguments leads to the following 
picture for the high-energy part of the muon flux: the 
production proceeds predominantly at stratospheric lev- 
els which are mainly subject to seasonal temperature 
(and thus air density) variations. The troposphere — with 
its daily modulation in temperature — only plays a minor 
role in determining the muon flux. The phenomenologi- 
cal relationship between muon flux and temperature can 
be written as 



AT, 



ax 



off 



(B-1) 



where (J^) is the (average) integral underground flux 
(with variation A/^) and ax is an effective tempera- 
ture coefficient with inferred values ar — 0.93 [20] and 
ar ^ 0.87 [25] for the LNGS and Soudan Mine sites, 
respectively. The effective temperature Tcs is obtained 
from the weighted average 



TcB — 



l^dXT{X)W{X) 
l^dXW{X) 



(B-2) 



where X is the atmospheric depth. The weight function 
W{X) depends on parameters such as the amount of in- 
clusive meson production in the forward fragmentation 
region, atmospheric attenuation parameters of mesons, 
and the muon spectral index. When computing Toff we 
use the expression for W given in [25] which includes the 
contribution of pions as well as kaons. 

The tight correlation between A/^ and AToff is firmly 
established by a whole range of independent measure- 
ments and at various geographic underground locations, 
such as at LNGS ^^M., at the Soudan Mine [21], or 
at the South Pole '[n]- For DAMA, the relevant LVD 
measurements are discussed in the introduction. Turn- 
ing to the CoGeNT site, the MINOS far detector has re- 
cently reported measurements on the underground muon 
flux spanning a total of five years [29]. Unfortunately, 
the published data does not overlap with the CoGeNT 
run starting Dec 4, 2009 and lasting 458 days. The re- 
sults are nevertheless relevant as they corroborate the 
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tight relation between effective temperature and mea- 
sured muon flux. In order to make statements about 
the muon background hypothesis for CoGeNT we have 
obtained cUmate data from the Integrated Global Ra- 
diosonde Archive (IGRA) for the nearby International 
Falls (MN) weather station [42]. We compute daily Tcs 
values by evaluating a discretized version of (B-2 1 includ- 



ing only high-quality data reaching atmospheric pressure 
levels below 40 hPa. The result of this procedure is shown 
in Fig. |12| which presents effective temperature values for 
a time span covering both, MINOS and CoGeNT, as in- 
dicated. We note in passing that the same IGRA data 
set was also used by the MINOS collaboration itself to 
confirm Teg values obtained from another climate archive 
with good agreement. 



Appendix C: Activation of long-lived isotopes 

A loop-hole that exists in the analysis of assessing the 
viability of the muon hypothesis for DAMA and CoGeNT 
is the following: What if one significantly breaks the 
temporal coincidence between muon arrival time in and 
around the detector and the time at which the signal 
is produced. "Significantly" here means a shift in time 
which at least on the order of weeks and more. It is clear 
then that one cannot compare the time series of muon 
flux data with the one from DAMA and CoGeNT at face 
value. 

The muon flux itself can give rise to cosmogenically ac- 
tivated isotopes. For example, the most inner shielding 
of the DAMA detector is made from massive copper so 
that for example a reaction ^^Cu -I- fi'^ — ^^Zn -I- i?^ with 
a delayed decay of ^^Zn with a half-life of 244 days is 
conceivable. Though the weak-scale cross section draws 
doubt on the efficiency of such a process, neutron capture 
reactions have cross sections many orders of magnitude 
larger. Indeed, it was one of the central points in Ref. [15] 
to point to the importance of (n, 7) reactions for back- 
ground considerations in direct detection experiments. 
Let us at this point not worry about the particular pro- 
cess but rather tackle the problem from the time series 
point of view. The question we are after is whether it 
is possible to induce a time-delay t in the signal such 
that it features a maximum in the year which is compat- 
ible with DAMA (and possibly by CoGenT.) As we have 



seen, the muon fiux peaks at least a month later than 
the signals under consideration. Therefore, one requires 
i > 300 days. 

For simplicity, let us model the muon flux as a 
sinosoidal function, /^(t) = /o[l -I- Im cos (wi)] where we 
are choosing the origin of time to coincide with its maxi- 
mum. Call n the number density of a stable element and 
n* its "excited" state formed upon a generalized interac- 
tion with the muon flux and which decays with a lifetime 
of r A~^. The abundance of n* is obtained by solving 



dn*{t) 
dt 



= -\n* (t) + 5o [1 + S,n cos {ujt)] (C-1) 
with an overall source term So with modulating fraction 
Sm for activation. In the physical picture, Sq = Xactn, 
where Aact is the rate for activation and given by the 
product of the activating (neutron) flux times the cross 
section (with a potential average over the flux energy 
spectrum.) Clearly, oc Iq holds with the same phase 
and modulation amplitude in Eq. (C-1) as in I^(t). 

If the events seen by DAMA (or CoGeNT) are some- 
how related to the decay of the radioactive isotope then 
the event rate in the detector will be proportional to 
R (X Xn*{t) which is just the usual relation for the num- 
ber of radioactive decays. The rate equation Eq. (C-1) 
can be solved exactly, 

n* (t) = Coe-^* + Sot + cos {cut - 0) , (C-2) 



where tan0 = uj/X and Sm — SoSm/ \/A^ + w^. The flrst 
term is transitory and vanishes for t ^ t. Therefore, the 
event rate will approach 



R-^ So + XSm cos[uj (t - t)]. 



(C-3) 



with the time delay, t = 6/uj. The latter is always a 
positive quantity with a minimum of zero days and max- 
imum of 91.25 days (quarter year). The phase delay is 
maximal {6 — w/2) when the lifetime of the decay is 
much longer than one year r ^ lyr/27r. Also, interest- 
ingly, the modulation amplitude is always smaller than 
the muon flux, since XS„i/So = (A/v'A^~-i-~w^)/m- Thus, 
we see that this muon related background is not optimal 
for explaining DAMA or CoGeNT, since the maximum 
of the signal can only be shifted into the fall, but not fur- 
ther into spring as would be required by DAMA and/or 
CoGeNT, and the modulation amplitude of the signal is 
always reduced from the muon's amplitude. 
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